Function-on-Scalar regressions for the white mold data

Author

Denis Shah

Published

March 4, 2025

Packages

library(fda)
library(refund)
library(fdatest)
library(tidyverse)
library(tictoc)
library(cowplot)
library(kableExtra)

Load and process the data

wm_load <- readr::read_csv(here::here("Data", "data_model_plus_weather_filtered.csv"), show_col_types = FALSE)

# Simplify things. Keep only the weather-related variables and others needed for calculations
wm_data <-
  wm_load %>%
  dplyr::select(subject, date, planting.date, sampling.date, wm, d2m, t2m_mean, t2m_max, t2m_min, st, sp, sm, rh) %>%
  # Do the filtering steps before doing any calculations or feature engineering:
  # wm_load has identical data for each of the sampling dates, which is why there is a filtering step on sampling.date.
  dplyr::group_by(subject) %>% 
  dplyr::filter(sampling.date == max(sampling.date)) %>%
  dplyr::ungroup() %>% 
  # Add the response variable (wm present or absent; binary):
  dplyr::group_by(subject) %>% 
  dplyr::mutate(wm = (mean(wm, na.rm = T) > 0)*1) %>% 
  # wm as a factor:
  dplyr::mutate(wm = factor(wm, levels = c(0, 1))) %>%
  dplyr::ungroup() %>% 
  dplyr::filter(!is.na(wm)) %>% 
  # Calculate dap (as a numeric):
  dplyr::mutate(dap = as.numeric(date - planting.date)) %>%
  # Convert temperatures from Kelvin to degrees Celsius:
  dplyr::mutate(across(d2m:st, ~ .x - 273.15)) %>%
  # dewpoint depression:
  dplyr::mutate(dpd = t2m_mean - d2m) %>%
  # surface pressure in kPa:
  dplyr::mutate(sp = sp/1000) %>%
  # Ratio of soil temperature to soil moisture:
  dplyr::mutate(stsm = st/sm) %>%
  # Difference in max and min temperatures:
  dplyr::mutate(tdiff = t2m_max - t2m_min) %>%
  # Calculate saturation vapor pressure (es):
  dplyr::mutate(es = 0.61078 * exp((17.27 * t2m_mean) / (t2m_mean + 237.3))) %>% 
  # Calculate actual vapor pressure (ea):
  dplyr::mutate(ea = (rh / 100) * es) %>% 
  # Calculate VPD (kPa):
  dplyr::mutate(vpd = es - ea) %>% 
  # estimate GDD (base 0). NB: base 0 is reasonable for snap bean; see the Jenni et al. (2000) paper. I think we want GDD to start accumulating the day after planting date onwards...
  # I use the day after planting, because we don't know exactly what time of the day the field was planted.
  dplyr::mutate(gddi = ifelse(dap <= 0, 0, (t2m_max + t2m_min)*0.5 - 0)) %>%
  dplyr::group_by(subject) %>% 
  # We don't want the gddi column after creating gdd:
  dplyr::mutate(gdd = cumsum(gddi), .keep = "unused") %>%
  dplyr::ungroup() %>%
  # Keep only the columns you need:
  dplyr::select(subject, dap, wm:gdd, -es, -ea) %>%
  # Filter the dataset to dap <= 50:
  dplyr::filter(dap <= 50) %>%
  # wm as a contrast: -1 = absent, 1 = present (needed for function-on-scalar regression):
  dplyr::mutate(wm = as.numeric(wm)-1) %>%
  dplyr::mutate(wm = 2*wm-1) 
data.prep <- function(x) {
  # Prepare the data for function-on-scalar regression
  # N x J matrix of the functional predictor 
  # Args:
  #  x = unquoted variable name
  # Returns:
  #  a list with the following named elements:
  #   y = N x J matrix of the functional response
  #   x = a vector of the wm status coded as -1, 1
  #   yind =  a vector of the evaluation days, which is of length 81 from -30 to 50
  
  .x <- enquo(x)
  
  # The environmental data:
  ez <-
    wm_data %>%
    dplyr::select(subject, dap, !!.x) %>%
    tidyr::pivot_wider(id_cols = subject, names_from = dap, names_prefix = "", values_from = !!.x) %>%
    dplyr::select(-subject) %>%
    as.matrix()
  
  colnames(ez) <- NULL
  
  # the vector of wm status:
  x.wm <- 
    wm_data %>%
    dplyr::group_by(subject) %>%
    dplyr::summarise(wm = mean(wm)) %>%
    dplyr::pull(wm)
  
  # the vector of the evaluation days:
  days <- seq(-30, 50)
  
  # the final data frame:
  df <- list(y = ez, x = x.wm, yind = days)
  
  return(df)
  
}

# Example of use:
# dat <- data.prep(d2m)

fANOVA <- function(dat) {
  # Performs a functional ANOVA and plots the estimated coefs for the overall mean and 
  # the difference between wm(0) and wm(1).
  #
  # Args:
  #   dat: prepared data from calling the data.prep function
  #
  # Returns:
  #   a ggplot graphic of the estimated beta(t) coefs.
  
  
  # k = 30 gives wiggliness, no oversmoothing
  m2 <- pffr(y~x, yind = dat$yind, data = dat,
             bs.yindex = list(bs = "ps", k = 30, m = c(2, 1)), 
             bs.int = list(bs = "ps", k = 30, m = c(2, 1)))
  
  
  # The smooth coefficients for the overall mean and beta(t):
  z <- coef(m2)$smterms
  
  # 1. For the overall mean:
  a1 <- z[[1]]$coef[, "yindex.vec"]
  a2 <- z[[1]]$coef[, "value"]
  a3 <- z[[1]]$coef[, "se"]
  
  # 2. For beta(t):
  b1 <- z[[2]]$coef[, "yindex.vec"]
  b2 <- z[[2]]$coef[, "value"]
  b3 <- z[[2]]$coef[, "se"]

  
  # Extract parts of the output to produce a more attractive plot in ggplot.
  # The coefs for the overall mean, adding an approx. 95% CI:
  X <- data.frame(x = a1, coef.mean = a2, lower.mean = a2 - 1.96*a3, upper.mean = a2 + 1.96*a3, coef = "Overall mean")
  
  # The coefs for the time trend of the mean difference between epidemics and non-epidemics:
  Y <- data.frame(x = b1, coef.mean = b2, lower.mean = b2 - 1.96*b3, upper.mean = b2 + 1.96*b3, coef = "Difference")
  
  Z <- rbind(X, Y)
  
  # To get different lines in each facet, you need another data.frame:
  hline.data <- data.frame(z = c(0), coef = c("Difference"))
  
  breaks <- seq(-30, 50, 10)
  labels <- seq(-30, 50, 10)
  
  Z %>%
  dplyr::mutate(coef = factor(coef, levels = c("Difference", "Overall mean"))) %>%
  ggplot(., aes(x = x, y = coef.mean)) +
  annotate("rect", ymin = -Inf, ymax = Inf, xmin = 35, xmax = 50, fill = "steelblue", alpha = 0.2) +
  geom_ribbon(aes(ymin = lower.mean, ymax = upper.mean), alpha = 0.2) +
  geom_path(linewidth = 1.2) +
  geom_hline(aes(yintercept = 0), color = "grey", linetype = "dashed") +
  geom_vline(aes(xintercept = 0), color = "grey", linetype = "dashed") +
  theme_bw() +
  facet_grid(coef ~ ., scales = "free_y") +
  theme(strip.text = element_text(face = "bold", size = rel(1.0))) +
  scale_x_continuous(name = "Days relative to sowing", breaks = breaks, labels = labels) +
  theme(axis.title.x = element_text(face = "bold", size = 11)) +
  ylab("Coefficient function") +
  theme(axis.title.y = element_text(face = "bold", size = 11))
  }

Function-on-scalar regression

Dew point

fANOVA(dat = data.prep(d2m))
using seWithMean for  s(yindex.vec) .

Mean air temperature

fANOVA(dat = data.prep(t2m_mean))
using seWithMean for  s(yindex.vec) .

Max air temperature

fANOVA(dat = data.prep(t2m_max))
using seWithMean for  s(yindex.vec) .

Min air temperature

fANOVA(dat = data.prep(t2m_min))
using seWithMean for  s(yindex.vec) .

Max - Min air temperature

fANOVA(dat = data.prep(tdiff))
using seWithMean for  s(yindex.vec) .

Temperature-Dewpoint depression

fANOVA(dat = data.prep(dpd))
using seWithMean for  s(yindex.vec) .

Soil temperature

fANOVA(dat = data.prep(st))
using seWithMean for  s(yindex.vec) .

Soil moisture

fANOVA(dat = data.prep(sm))
using seWithMean for  s(yindex.vec) .

soil temperature:soil moisture ratio

fANOVA(dat = data.prep(stsm))
using seWithMean for  s(yindex.vec) .

Surface pressure

fANOVA(dat = data.prep(sp))
using seWithMean for  s(yindex.vec) .

Relative humidity

fANOVA(dat = data.prep(rh))
using seWithMean for  s(yindex.vec) .

VPD

fANOVA(dat = data.prep(vpd))
using seWithMean for  s(yindex.vec) .

Growing degree days

fANOVA(dat = data.prep(gdd))
using seWithMean for  s(yindex.vec) .


# We will illustrate fdatest with the d2m variable.

# For fdatest, we need to create two separate matrices (for wm = 0 and wm = 1), which we can then put into a list for convenience.
wm0 <- 
  wm_data %>%
  dplyr::filter(wm == -1) %>%
  dplyr::select(subject, dap, d2m) %>%
  tidyr::pivot_wider(id_cols = subject, names_from = dap, names_prefix = "", values_from = d2m) %>%
  dplyr::select(-subject) %>%
  as.matrix()
# and ...
colnames(wm0) <- NULL

wm1 <- 
  wm_data %>%
  dplyr::filter(wm == 1) %>%
  dplyr::select(subject, dap, d2m) %>%
  tidyr::pivot_wider(id_cols = subject, names_from = dap, names_prefix = "", values_from = d2m) %>%
  dplyr::select(-subject) %>%
  as.matrix()
colnames(wm1) <- NULL

# Performing the ITP:
ITP.result <- fdatest::ITP2bspline(wm0, wm1, B = 100)

# The function generates a print line for each iteration. To suppress that, wrap within sink:
tic()
{ sink(type = "message"); ITP.result <- fdatest::ITP2bspline(wm0, wm1, B = 1000); sink() }
toc()  # 9.18 sec

# Plotting the results of the ITP:
# (there are two plots. The first is of the individual curves. The 2nd is of the adjusted p-values)
plot(ITP.result, main = NULL, xrange = c(1, 365), xlab = 'Day')

# Plotting the p-values heatmap
# (I'm not finding it all that telling)
ITPimage(ITP.result, abscissa.range = c(0, 1))

# Selecting the significant components at 5% level:
which(ITP.result$corrected.pval < 0.05)

# Which corresponds to the following days (relative to sowing, where sowing is day = 0):
seq(-30, 50)[31:32]

# NEXT:
# Take the above code, and wrap into functions to (i) prep the data for input into fdatest, (ii) calling the fdatest functions to estimate the adjusted p values
make_kable <- function(...) {
  # kable and kableExtra styling to avoid repetitively calling the styling over and over again
  # See: https://stackoverflow.com/questions/73718600/option-to-specify-default-kableextra-styling-in-rmarkdown
  # knitr::kable(...) %>%
  kable(..., format = "html", row.names = TRUE, align = 'l') %>%
    kable_styling(bootstrap_options = c("striped"), position = "left", font_size = 11, full_width = FALSE) 
}

do.fdatest <- function(x) {
  # Performs an Interval Testing Procedure for testing the difference between the two functional wm populations evaluated on a uniform grid
  #
  # Args:
  #  x = unquoted variable name, e.g. d2m
  #
  # Returns:
  #   a tibble of the variable with a list vector of the days (relative to sowing) where the two populations differ functionally
  
  # Data prep for input to fdatest
  .x <- enquo(x)
  # wm0 = N x J matrix of the functional data for wm absent
  wm0 <- 
    wm_data %>%
    dplyr::filter(wm == -1) %>%
    dplyr::select(subject, dap, !!.x) %>%
    tidyr::pivot_wider(id_cols = subject, names_from = dap, names_prefix = "", values_from = !!.x) %>%
    dplyr::select(-subject) %>%
    as.matrix()
    # and ...
    colnames(wm0) <- NULL
    
  # wm1 = N x J matrix of the functional data for wm present  
  wm1 <- 
    wm_data %>%
    dplyr::filter(wm == 1) %>%
    dplyr::select(subject, dap, !!.x) %>%
    tidyr::pivot_wider(id_cols = subject, names_from = dap, names_prefix = "", values_from = !!.x) %>%
    dplyr::select(-subject) %>%
    as.matrix()
    colnames(wm1) <- NULL
  
  dat <- list(wm0 = wm0, wm1 = wm1)
  
  # The function generates a print line for each iteration. To suppress that, wrap within sink:
  # (set a seed for reproducibility)
  { sink(nullfile()); set.seed(86754309); foo <- ITP2bspline(dat$wm0, dat$wm1, B = 10000); sink() }

  # Selecting the significant components at 5% level:
  # Which corresponds to the following days (relative to sowing, where sowing is day = 0):
  z <- seq(-30, 50)[which(foo$corrected.pval < 0.05)]
  
  res <- tibble(var = rlang::as_name(.x), days = list(z))
  
  return(res)
}

# Example of use:
# do.fdatest(x = d2m)

filter.iwt <- function(x) {
  # Filter the interval-wise testing results to see the days that were significant
  # Args:
  #  x = the series (quoted character string)
  # Returns:
  #  a vector of the days where the series was different between wm = 0 and wm = 1
  iwt |>
  dplyr::filter(var == x) |>
  purrr::pluck("days", 1)
}


zee <- function(series) {
  # Output the start and end days of significant windows within a time series
  # Args:
  #  series = quoted character name of the series
  # Returns:
  #  a table
  v <- 
    filter.iwt(x = series) %>% 
    split(., cumsum(c(1, diff(.) != 1)))
  
  # Create an empty data frame with three named columns:
  z <- data.frame(matrix(
    vector(), 0, 3, dimnames = list(c(), c("series", "start", "end"))), 
    stringsAsFactors = F)
  
  # Now loop over v to pick out the start and end of the continuous windows:
  for (i in 1:length(v)) {
    b <- purrr::pluck(v, i) %>%  dplyr::first()
    c <- purrr::pluck(v, i) %>%  dplyr::last()
    z[i, "series"] <- series
    z[i, "start"] <- b
    z[i, "end"] <- c
    } # end for loop
  return(z)
}

# Examples of use:
# zee("t2m_mean")
# zee("d2m")

#  Example of the workflow:
###---###
# tst.d2m <- do.fdatest(x = d2m)

# You could output this way:
# pluck(iwt.list, "tst.d2m") %>% 
#   make_kable()

# But this shows the windows in a cleaner format:  
# zee("d2m") %>% 
  # make_kable()
###---###

# Now we're ready to roll...
load(here::here("FunctionalDataAnalysis", "FunctiononScalar", "iwt.RData"))

Interval-wise tests

Perform the interval-wise tests

tst.d2m <- do.fdatest(x = d2m)

tst.t2m_mean <- do.fdatest(x = t2m_mean)

tst.t2m_max <- do.fdatest(x = t2m_max)

tst.t2m_min <- do.fdatest(x = t2m_min)

tst.tdiff <- do.fdatest(x = tdiff)

tst.dpd <- do.fdatest(x = dpd)

tst.st <- do.fdatest(x = st)

tst.sm <- do.fdatest(x = sm)

tst.stsm <- do.fdatest(x = stsm)

tst.sp <- do.fdatest(x = sp)

tst.rh <- do.fdatest(x = rh)

tst.vpd <- do.fdatest(x = vpd)

tst.gdd <- do.fdatest(x = gdd)

I’ve already done the interval-wise tests, so the chunks in this section will just output the days that the test found to be significant.

Dew point

zee("d2m") %>% 
  make_kable()
series start end
1 d2m 0 0

Mean air temperature

zee("t2m_mean") %>% 
  make_kable()
series start end
1 t2m_mean 0 4

Max air temperature

zee("t2m_max") %>% 
  make_kable()
series start end
1 t2m_max 0 4

Min air temperature

zee("t2m_min") %>% 
  make_kable()
series start end
1 t2m_min 0 1

Max - Min air temperature

zee("tdiff") %>% 
  make_kable()
series start end
1 tdiff NA NA

Temperature-Dewpoint depression

zee("dpd") %>% 
  make_kable()
series start end
1 dpd NA NA

Soil temperature

zee("st") %>% 
  make_kable()
series start end
1 st 1 3

Soil moisture

zee("sm") %>% 
  make_kable()
series start end
1 sm -4 3
2 sm 5 15
3 sm 17 24
4 sm 40 49

soil temperature:soil moisture ratio

zee("stsm") %>% 
  make_kable()
series start end
1 stsm 24 25
2 stsm 35 44

Surface pressure

zee("sp") %>% 
  make_kable()
series start end
1 sp NA NA

Relative humidity

zee("rh") %>% 
  make_kable()
series start end
1 rh NA NA

VPD

zee("vpd") %>% 
  make_kable()
series start end
1 vpd NA NA
# I want to put all the data frames into a list (to pass to bind_rows). Found the solution at:
# https://stackoverflow.com/questions/26738302/make-list-of-objects-in-global-environment-matching-certain-string-pattern
# so don't have to write them out one-by-one
Pattern1 <- grep("^tst.", names(.GlobalEnv), value = TRUE)
iwt.list <- do.call("list", mget(Pattern1))
# Don't need the tst. objects no more:
# rm(list = Pattern1)

# Now bind all the data frames into one:
iwt <- dplyr::bind_rows(iwt.list)

# Save iwt: 
save(iwt.list, iwt, file = here::here("FunctionalDataAnalysis", "FunctiononScalar", "iwt.RData"))

Summary of the test results

# Let's use the filter.iwt function to have a look at the intervals that were significant (all series, ignoring gdd):
# air temperature-related
filter.iwt(x = "t2m_mean")
filter.iwt(x = "t2m_max")
filter.iwt(x = "t2m_min")
filter.iwt(x = "tdiff")
filter.iwt(x = "d2m") 
filter.iwt(x = "dpd")

# soil-related
filter.iwt(x = "sm")
filter.iwt(x = "st")
filter.iwt(x = "stsm")

# surface pressure
filter.iwt(x = "sp")

# moisture-related
filter.iwt(x = "rh")
filter.iwt(x = "vpd")

# We want to be able to cut the significant periods into continuous windows should those exist.
# The code to do so is adapted from here:
# https://stackoverflow.com/questions/5222061/create-grouping-variable-for-consecutive-sequences-and-split-vector
# Examples:
filter.iwt(x = "t2m_max") %>% split(., cumsum(c(1, diff(.) != 1)))
filter.iwt(x = "stsm") %>% split(., cumsum(c(1, diff(.) != 1)))


# What we're ultimately interested in is parsing the significant dates into any continuous segments.
# The find the start and end of these segments, so that we can input those to a dataframe for plotting. For that we use the function `zee` (see the functions chunk)

# Create a vector of all the series names:
e <- c("t2m_mean", "t2m_max", "t2m_min", "tdiff", "d2m", "dpd", "sm", "st", "stsm", "sp", "rh", "vpd")

# May over e and bind the rows:
g <- purrr::map(e, zee) |> list_rbind()
# Code for plotting the segments derived from here:
# https://stackoverflow.com/questions/35322919/grouped-data-by-factor-with-geom-segment

# Assign variable labels:
var_labels <- tribble(
  ~variable,      ~variable_label, 
  "d2m",       "Dewpoint (°C)",
  "t2m_mean",  "Mean air temperature (°C)",
  "t2m_max",   "Max. air temperature (°C)",
  "t2m_min",   "Min. air temperature (°C)",
  "tdiff",     "Max. - Min. air temperature (°C)",
  "dpd",       "Dew point depression (°C)",
  "st",        "Soil temperature (°C)",
  "sm",        "Soil moisture (m³/m³)",
  "stsm",      "Soil temperature:moisture ratio",
  "sp",        "Surface pressure (kPa)",
  "rh",        "Relative humidity (%)", 
  "vpd",       "Vapor pressure deficit (kPa)"
) %>% 
  tibble::deframe()

g %>% 
  ggplot(., aes(ymin = start, ymax = end, x = series)) + 
  # Changing the order on the x-axis for the categories:
  scale_x_discrete(limits = rev(c("t2m_mean", "t2m_max", "t2m_min", "tdiff", "d2m", "dpd", "sm", "st", "stsm", "sp", "rh", "vpd")), labels = var_labels, name = NULL) + 
  # Pay attention to layers. Do coord_flip first, then annotations to appear in the background, segment the topmost layer:
  coord_flip() +
  # Add some guides:
  geom_hline(aes(yintercept = 0), color = "grey", linetype = "dashed") +
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = 35, ymax = 50, fill = "steelblue", alpha = 0.2) +
    geom_linerange(colour = "grey20", position = position_dodge(width = 0.2), linewidth = 3, na.rm = TRUE) + 
    # theme_bw() +
    theme_half_open(font_size = 12)+
    labs(x = "Series", y = "Days relative to sowing") +
    theme(axis.title = element_text(face = "bold", size = 11))

ggsave("figs/fda_difference_period.png", dpi = 600, bg = "white")
weather_vars <- function(whichvar, start, end) {
  # Calculate the weather-based summary variable for each observation
  # Args:
  #  whichvar = unquoted character string of the variable (series), e.g., sm
  #  start = the start day relative to sowing
  #  end = the end day of the window relative to sowing
  # Returns:
  #  a data frame with two columns (subject, the weather-based summary variable)
  
  var <- enquo(whichvar)
  # Create a name for the column to hold the summary variable:
  var_name <- paste(rlang::as_name(var), start, end, sep = "_")
  
  wm_data %>%
    dplyr::select(subject, dap, !!var) %>%
    dplyr::filter(dap %in% c(start:end)) %>%
    dplyr::group_by(subject) %>%
    dplyr::summarise("{var_name}" := mean(!!var))
  } # end of function


# These are the variables and windows we'll create summaries for:
# t2m_mean  start = 0, end = 4
# sm start = -4, end = 3
# sm start = 5, end = 15
# sm start = 17, end = 24
# sm start = 40, end = 49
# stsm start = 35, end = 44

u <- weather_vars(whichvar = t2m_mean, start = 0, end = 4)
v <- weather_vars(whichvar = sm, start = -4, end = 3)
w <- weather_vars(whichvar = sm, start = 5, end = 15)
x <- weather_vars(whichvar = sm, start = 17, end = 24)
y <- weather_vars(whichvar = sm, start = 40, end = 49)
z <- weather_vars(whichvar = stsm, start = 35, end = 44)

weather.vars <- purrr::reduce(list(u, v, w, x, y, z), dplyr::left_join, by = "subject")

# Save the weather.vars data frame:
save(weather.vars, file = "WeatherVars.RData")

Session Info

sessionInfo()
R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)

Matrix products: default


locale:
[1] LC_COLLATE=Portuguese_Brazil.utf8  LC_CTYPE=Portuguese_Brazil.utf8   
[3] LC_MONETARY=Portuguese_Brazil.utf8 LC_NUMERIC=C                      
[5] LC_TIME=Portuguese_Brazil.utf8    

time zone: America/Sao_Paulo
tzcode source: internal

attached base packages:
[1] splines   stats     graphics  grDevices datasets  utils     methods  
[8] base     

other attached packages:
 [1] kableExtra_1.4.0 cowplot_1.1.3    tictoc_1.2.1     lubridate_1.9.3 
 [5] forcats_1.0.0    stringr_1.5.1    dplyr_1.1.4      purrr_1.0.2     
 [9] readr_2.1.5      tidyr_1.3.1      tibble_3.2.1     ggplot2_3.5.1   
[13] tidyverse_2.0.0  fdatest_2.1.1    refund_0.1-37    fda_6.2.0       
[17] deSolve_1.40     fds_1.8          RCurl_1.98-1.16  rainbow_3.8     
[21] pcaPP_2.0-5      MASS_7.3-60.2    knitr_1.48      

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1   hdrcde_3.4         viridisLite_0.4.2  RLRsim_3.1-8      
 [5] farver_2.1.2       bitops_1.0-9       fastmap_1.2.0      pracma_2.4.4      
 [9] digest_0.6.37      timechange_0.3.0   lifecycle_1.0.4    cluster_2.1.6     
[13] magrittr_2.0.3     compiler_4.4.1     rlang_1.1.4        tools_4.4.1       
[17] utf8_1.2.4         yaml_2.3.10        labeling_0.4.3     grpreg_3.5.0      
[21] htmlwidgets_1.6.4  bit_4.5.0          mclust_6.1.1       here_1.0.1        
[25] xml2_1.3.6         abind_1.4-8        KernSmooth_2.23-24 gamm4_0.2-6       
[29] withr_3.0.2        grid_4.4.1         fansi_1.0.6        colorspace_2.1-1  
[33] scales_1.3.0       cli_3.6.3          mvtnorm_1.3-2      crayon_1.5.3      
[37] rmarkdown_2.28     ragg_1.3.3         generics_0.1.3     rstudioapi_0.17.0 
[41] tzdb_0.4.0         magic_1.6-1        minqa_1.2.8        parallel_4.4.1    
[45] vctrs_0.6.5        boot_1.3-30        Matrix_1.7-0       jsonlite_1.8.9    
[49] hms_1.1.3          bit64_4.5.2        systemfonts_1.1.0  glue_1.8.0        
[53] nloptr_2.1.1       stringi_1.8.4      gtable_0.3.5       lme4_1.1-35.5     
[57] munsell_0.5.1      pillar_1.9.0       htmltools_0.5.8.1  R6_2.5.1          
[61] textshaping_0.4.0  ks_1.14.3          rprojroot_2.0.4    vroom_1.6.5       
[65] evaluate_1.0.1     lattice_0.22-6     highr_0.11         renv_1.1.2        
[69] pbs_1.1            Rcpp_1.0.13        svglite_2.1.3      nlme_3.1-164      
[73] mgcv_1.9-1         xfun_0.48          pkgconfig_2.0.3